LOAD DATA

df <- readRDS("../data/models/social-risk-crash-rate-data.rds")

HOT-ENCODE YEAR

df$year <- as.factor(df$year)
year_dummies <- model.matrix(~ year - 1, data = df)

df <- cbind(df[ , !(names(df) %in% c("year"))], year_dummies)

SELECT TRANSFORMATION FOR TARGET

compare_transformations <- function(df, target_var, features) {
  
  # Load required packages
  if (!require("MASS")) install.packages("MASS")
  if (!require("bestNormalize")) install.packages("bestNormalize")
  library(MASS)
  library(bestNormalize)
  
  y_raw <- df[[target_var]]
  X <- df[, features, drop = FALSE]
  X[] <- lapply(X, function(col) as.numeric(as.character(col)))
  
  # Log Transformation 
  y_log <- log1p(y_raw)
  lm_log <- lm(y_log ~ ., data = data.frame(y_log, X))
  pred_log <- predict(lm_log, newdata = data.frame(X))
  rmse_log <- sqrt(mean((y_log - pred_log)^2))
  
  # Box-Cox Transformation 
  y_bc_input <- ifelse(y_raw <= 0, 0.000001, y_raw)
  
  bc <- MASS::boxcox(y_bc_input ~ ., data = data.frame(y_bc_input, X), lambda = seq(-2, 2, 0.1))
  best_lambda <- bc$x[which.max(bc$y)]
  y_bc <- if(best_lambda == 0) {
    log(y_bc_input)
  } else {
    (y_bc_input^best_lambda - 1)/best_lambda
  }
  
  lm_bc <- lm(y_bc ~ ., data = data.frame(y_bc, X))
  pred_bc <- predict(lm_bc, newdata = data.frame(X))
  rmse_bc <- sqrt(mean((y_bc - pred_bc)^2))
  
  # Yeo-Johnson Transformation
  yj_obj <- yeojohnson(y_raw)
  y_yj <- predict(yj_obj)
  
  lm_yj <- lm(y_yj ~ ., data = data.frame(y_yj, X))
  pred_yj <- predict(lm_yj, newdata = data.frame(X))
  rmse_yj <- sqrt(mean((y_yj - pred_yj)^2))
  
  results <- data.frame(
    Transformation = c("Log (log1p)", "Box-Cox", "Yeo-Johnson"),
    RMSE = c(rmse_log, rmse_bc, rmse_yj),
    Best_Lambda_BoxCox = c(NA, best_lambda, NA)
  )
  
  print(results)
  
  return(list(
    log_model = lm_log,
    boxcox_model = lm_bc,
    yj_model = lm_yj,
    boxcox_lambda = best_lambda,
    transformation_results = results
  ))
}
# Select target variable
target_var <- "crash_rate_per_1000"

df <- df %>% dplyr::select(-injury_rate_per_1000, -fatality_rate_per_1000,  -total_population, -geoid)

features <- setdiff(names(df), c("crash_rate_per_1000", "borough"))
comparison <- compare_transformations(df, "crash_rate_per_1000", features)
NA

LOG TRANSFORM TARGET

X <- df %>% dplyr::select(-target_var, -borough)

# Log-transform the target variable
y <- log1p(df[[target_var]])  

hist(y, breaks = 50, main = "Log-Transformed Crash Rate", xlab = "log1p(crash_rate_per_1000)")

HYPERPARAMETER TUNING with Optuna

## CONVERT TO DMATRIX
dtrain_all <- xgb.DMatrix(data = as.matrix(X), label = y)

# Define Spatial Folds 
boroughs <- unique(df$borough)
folds <- lapply(boroughs, function(b) which(df$borough != b))
names(folds) <- boroughs

## Start Python venv
reticulate::use_virtualenv("r-reticulate", required = TRUE)

## OPTUNA-BASED SPATIAL CV
optuna <- import("optuna")

# Objective Function 
objective <- function(trial) {
  params <- list(
    booster = "gbtree",
    objective = "reg:squarederror",
    eval_metric = "rmse",
    eta = trial$suggest_float("eta", 0.01, 0.3, log = TRUE),
    max_depth = trial$suggest_int("max_depth", 3, 12),
    min_child_weight = trial$suggest_int("min_child_weight", 1, 10),
    subsample = trial$suggest_float("subsample", 0.5, 1.0),
    colsample_bytree = trial$suggest_float("colsample_bytree", 0.5, 1.0),
    gamma = trial$suggest_float("gamma", 0, 10),
    lambda = trial$suggest_float("lambda", 0, 10),
    alpha = trial$suggest_float("alpha", 0, 10)
  )
  
  cv_results <- xgb.cv(
    params = params,
    data = dtrain_all,
    nrounds = 150,
    folds = folds,
    early_stopping_rounds = 20,
    verbose = 0,
    nthread = detectCores() - 1
  )
  
  return(cv_results$evaluation_log$test_rmse_mean[cv_results$best_iteration])
}


set.seed(2025)
study <- optuna$create_study(direction = "minimize")
study$optimize(objective, n_trials = 50)

best_params <- study$best_params
best_params$objective <- "reg:squarederror"
print(best_params)

$eta [1] 0.2586295

$max_depth [1] 12

$min_child_weight [1] 4

$subsample [1] 0.9921435

$colsample_bytree [1] 0.5823858

$gamma [1] 0.08316905

$lambda [1] 4.260907

$alpha [1] 0.2049112

TRAIN/TEST SPLIT

# Set seed
set.seed(2025)

# Split by index
train_index <- createDataPartition(y, p = 0.8, list = FALSE)
X_train <- X[train_index, ]
y_train <- y[train_index]
X_test <- X[-train_index, ]
y_test <- y[-train_index]

# Convert to xgb.DMatrix
dtrain <- xgb.DMatrix(data = as.matrix(X_train), label = y_train)
dtest <- xgb.DMatrix(data = as.matrix(X_test), label = y_test)

TRAIN MODEL

# Set seed
set.seed(2025)

# Training with parallel processing
final_model <- xgb.train(
  params = list(
    eta = best_params$eta,
    max_depth = best_params$max_depth,
    gamma = best_params$gamma,
    colsample_bytree = best_params$colsample_bytree,
    min_child_weight = best_params$min_child_weight,
    subsample = best_params$subsample,
    objective = "reg:squarederror",
    eval_metric = "rmse"
  ),
  data = dtrain,
  nrounds = 1000,
  watchlist = list(train = dtrain, test = dtest),
  early_stopping_rounds = 20,
  verbose = 1,
  nthread = detectCores() - 1
)

SAVE MODEL

# Create directory if it doesn't exist
if (!dir.exists("../data/models")) {
  dir.create("../data/models", recursive = TRUE)
}

# Save the final XGBoost model
saveRDS(final_model, file = "../data/models/xgb_model.rds")

# Save the best parameters
saveRDS(best_params, file = "../data/models/xgb_best_params.rds")

cat("Model and parameters saved to ../data/models/")
# Load the final XGBoost model
final_model <- readRDS("../data/models/xgb_model.rds")

# Load the best parameters
best_params <- readRDS("../data/models/xgb_best_params.rds")

MODEL EVALUATION

set.seed(2025)

# Predict on test set
preds <- predict(final_model, as.matrix(X_test))

# Metrics 
rmse <- sqrt(mean((y_test - preds)^2))
mae <- mean(abs(y_test - preds))
r2 <- 1 - (sum((y_test - preds)^2) / sum((y_test - mean(y_test))^2))

cat("Model Evaluation Metrics:\n")
Model Evaluation Metrics:
cat("  RMSE:", rmse, "\n")
  RMSE: 0.3396146 
cat("  MAE :", mae, "\n")
  MAE : 0.2422538 
cat("  R²  :", r2, "\n\n")
  R²  : 0.425285 
# Residuals 
residuals <- y_test - preds
residual_df <- data.frame(
  actual = y_test,
  predicted = preds,
  residuals = residuals
)

# Predicted vs Actual
p1 <- residual_df %>%
  ggplot(aes(x = actual, y = predicted)) +
  geom_point(alpha = 0.5) +
  geom_abline(slope = 1, intercept = 0, color = "red") +
  theme_minimal() +
  labs(title = "Predicted vs Actual Crash Rates",
       x = "Actual",
       y = "Predicted")

# Residuals vs Predicted 
p2 <- residual_df %>%
  ggplot(aes(x = predicted, y = residuals)) +
  geom_point(alpha = 0.5, color = "blue") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  theme_minimal() +
  labs(title = "Residuals vs Predicted",
       x = "Predicted",
       y = "Residuals")

# Residual Density 
p3 <- ggplot(residual_df, aes(x = residuals)) +
    geom_histogram(
        aes(y = ..density..),  
        binwidth = 0.05, 
        fill = "steelblue", 
        color = "white",
        alpha = 0.7
    ) +
    geom_density(
        color = "red", 
        linewidth = 1, 
        adjust = 1.5  
    ) +
    theme_minimal() +
    labs(
        title = "Residual Distribution", 
        x = "Residuals", 
        y = "Density"   
    )


ggsave("../report/plots/predicted_vs_actual_values_plot.png", p1, width = 10, height = 6, dpi = 300)
ggsave("../report/plots/resisuals_vs_predicted_values_plot.png", p2, width = 10, height = 6, dpi = 300)
ggsave("../report/plots/residual_density_plot.png", p3, width = 10, height = 6, dpi = 300)

# Grid Plots
grid_plot <- p1 + p2 + p3 + plot_layout(ncol = 3)

ggsave(
  "../report/plots/residuals_diagnostic_grid.png",
  grid_plot,
  width = 18, 
  height = 6,
  dpi = 300
)

Model Evaluation Metrics: RMSE: 0.3396146 MAE : 0.2422538 R² : 0.425285 ## SHAP EXPLAINABILITY

# Compute SHAP values
shap_values <- shap.values(xgb_model = final_model, X_train = as.matrix(X_train))
shap_long <- shap.prep(shap_contrib = shap_values$shap_score, X_train = as.matrix(X_train))

# SHAP summary plot
print(shap.plot.summary(shap_long))

shap <- shap.plot.summary(shap_long)
ggsave("../report/plots/shap_summary_plot.png", shap, width = 10, height = 6, dpi = 300)

INSPECT TREE ENSEMBLE

xgb.plot.multi.trees(model = final_model)
# SHAP Dependence and Interaction Plots

# Interaction values
shap_interaction_values <- predict(
  final_model,
  as.matrix(X_train),
  predinteraction = TRUE
)

saveRDS(shap_interaction_values, "../data/models/shap_interactions.rds")
# Define the specific top features we want to analyze
selected_features <- c(
  "post_pandemic",
  "pct_age_65_plus",
  "pct_graduate_degree",
  "pct_in_labor_force",
  "median_income",
  "pct_commute_medium"
)

# Create pretty names for the selected features
feature_names <- c(
  "Pre vs Post Pandemic",
  "% Population 65+ Years",
  "% with Graduate Degree",
  "% in Labor Force",
  "Median Income",
  "% with Medium Commute"
)

# Verify all selected features exist in X_train
missing_features <- setdiff(selected_features, colnames(X_train))
if (length(missing_features) > 0) {
  stop(paste("The following features are missing from X_train:", 
             paste(missing_features, collapse = ", ")))
}

# Generate PDP plots with error handling
pdp_plots <- map2(selected_features, feature_names, function(f, name) {
  tryCatch({
    pd <- pdp::partial(
      final_model, 
      pred.var = f, 
      train = as.matrix(X_train), 
      grid.resolution = 30
    )
    
    autoplot(pd) +
      ggtitle(name) +
      theme_minimal() +
      theme(
        axis.title = element_blank(),
        axis.text = element_text(size = 8),
        plot.title = element_text(size = 10, hjust = 0.5),
        plot.margin = unit(c(1, 1, 1, 1), "lines")
      )
  }, error = function(e) {
    message(paste("Error generating PDP for", f, ":", e$message))
    return(NULL)
  })
})

# Remove any NULL entries from failed plots
pdp_plots <- purrr::compact(pdp_plots)

# Create directory if it doesn't exist
if (!dir.exists("../report/plots/pdp")) {
  dir.create("../report/plots/pdp", recursive = TRUE)
}

# Save individual plots
walk2(pdp_plots, selected_features, ~{
  if (!is.null(.x)) {
    ggsave(
      filename = paste0("../report/plots/pdp/pdp_", .y, ".png"),
      plot = .x,
      width = 5,
      height = 4,
      dpi = 300
    )
  }
})

# Create and save grid plot
if (length(pdp_plots) > 0) {
  grid_plot <- gridExtra::grid.arrange(
    grobs = pdp_plots, 
    ncol = 2,
    top = "Partial Dependence Plots for Key Features"
  )
  
  ggsave(
    filename = "../report/plots/pdp/pdp_grid.png",
    plot = grid_plot,
    width = 10,
    height = 12,
    dpi = 300
  )
  
  # Optionally display the grid plot
  print(grid_plot)
} else {
  warning("No PDP plots were successfully generated.")
}
LS0tCnRpdGxlOiAiWEdCb29zdCBNb2RlbGluZyBSIE5vdGVib29rIgphdXRob3I6ICJBSiBTdHJhdW1hbi1TY290dCIKZGF0ZTogImByIFN5cy5EYXRlKClgIgpvdXRwdXQ6CiAgaHRtbF9kb2N1bWVudDoKICAgIGRmX3ByaW50OiBwYWdlZAogIGh0bWxfbm90ZWJvb2s6IGRlZmF1bHQKICBwZGZfZG9jdW1lbnQ6IGRlZmF1bHQKLS0tCgpgYGB7ciBzZXR1cCwgaW5jbHVkZT1GQUxTRX0Ka25pdHI6Om9wdHNfY2h1bmskc2V0KGVjaG8gPSBUUlVFLCB3YXJuaW5nID0gRkFMU0UsIG1lc3NhZ2UgPSBGQUxTRSkKbGlicmFyeSh0aWR5dmVyc2UpCmxpYnJhcnkoY2FyZXQpCmxpYnJhcnkoeGdib29zdCkKbGlicmFyeShTSEFQZm9yeGdib29zdCkKbGlicmFyeShNYXRyaXgpCmxpYnJhcnkocmV0aWN1bGF0ZSkKbGlicmFyeShkb1BhcmFsbGVsKQpsaWJyYXJ5KE1ldHJpY3MpCmxpYnJhcnkocGRwKQpsaWJyYXJ5KERBTEVYKQpsaWJyYXJ5KGdyaWRFeHRyYSkKbGlicmFyeShwYXRjaHdvcmspCmBgYAoKIyMgTE9BRCBEQVRBCgpgYGB7ciBsb2FkLWRhdGF9CmRmIDwtIHJlYWRSRFMoIi4uL2RhdGEvbW9kZWxzL3NvY2lhbC1yaXNrLWNyYXNoLXJhdGUtZGF0YS5yZHMiKQpgYGAKCiMjIEhPVC1FTkNPREUgYFlFQVJgCgpgYGB7ciBob3QtZW5jb2RlLXllYXJ9CmRmJHllYXIgPC0gYXMuZmFjdG9yKGRmJHllYXIpCnllYXJfZHVtbWllcyA8LSBtb2RlbC5tYXRyaXgofiB5ZWFyIC0gMSwgZGF0YSA9IGRmKQoKZGYgPC0gY2JpbmQoZGZbICwgIShuYW1lcyhkZikgJWluJSBjKCJ5ZWFyIikpXSwgeWVhcl9kdW1taWVzKQpgYGAKCiMjIFNFTEVDVCBUUkFOU0ZPUk1BVElPTiBGT1IgVEFSR0VUCgpgYGB7cn0KY29tcGFyZV90cmFuc2Zvcm1hdGlvbnMgPC0gZnVuY3Rpb24oZGYsIHRhcmdldF92YXIsIGZlYXR1cmVzKSB7CiAgCiAgIyBMb2FkIHJlcXVpcmVkIHBhY2thZ2VzCiAgaWYgKCFyZXF1aXJlKCJNQVNTIikpIGluc3RhbGwucGFja2FnZXMoIk1BU1MiKQogIGlmICghcmVxdWlyZSgiYmVzdE5vcm1hbGl6ZSIpKSBpbnN0YWxsLnBhY2thZ2VzKCJiZXN0Tm9ybWFsaXplIikKICBsaWJyYXJ5KE1BU1MpCiAgbGlicmFyeShiZXN0Tm9ybWFsaXplKQogIAogIHlfcmF3IDwtIGRmW1t0YXJnZXRfdmFyXV0KICBYIDwtIGRmWywgZmVhdHVyZXMsIGRyb3AgPSBGQUxTRV0KICBYW10gPC0gbGFwcGx5KFgsIGZ1bmN0aW9uKGNvbCkgYXMubnVtZXJpYyhhcy5jaGFyYWN0ZXIoY29sKSkpCiAgCiAgIyBMb2cgVHJhbnNmb3JtYXRpb24gCiAgeV9sb2cgPC0gbG9nMXAoeV9yYXcpCiAgbG1fbG9nIDwtIGxtKHlfbG9nIH4gLiwgZGF0YSA9IGRhdGEuZnJhbWUoeV9sb2csIFgpKQogIHByZWRfbG9nIDwtIHByZWRpY3QobG1fbG9nLCBuZXdkYXRhID0gZGF0YS5mcmFtZShYKSkKICBybXNlX2xvZyA8LSBzcXJ0KG1lYW4oKHlfbG9nIC0gcHJlZF9sb2cpXjIpKQogIAogICMgQm94LUNveCBUcmFuc2Zvcm1hdGlvbiAKICB5X2JjX2lucHV0IDwtIGlmZWxzZSh5X3JhdyA8PSAwLCAwLjAwMDAwMSwgeV9yYXcpCiAgCiAgYmMgPC0gTUFTUzo6Ym94Y294KHlfYmNfaW5wdXQgfiAuLCBkYXRhID0gZGF0YS5mcmFtZSh5X2JjX2lucHV0LCBYKSwgbGFtYmRhID0gc2VxKC0yLCAyLCAwLjEpKQogIGJlc3RfbGFtYmRhIDwtIGJjJHhbd2hpY2gubWF4KGJjJHkpXQogIHlfYmMgPC0gaWYoYmVzdF9sYW1iZGEgPT0gMCkgewogICAgbG9nKHlfYmNfaW5wdXQpCiAgfSBlbHNlIHsKICAgICh5X2JjX2lucHV0XmJlc3RfbGFtYmRhIC0gMSkvYmVzdF9sYW1iZGEKICB9CiAgCiAgbG1fYmMgPC0gbG0oeV9iYyB+IC4sIGRhdGEgPSBkYXRhLmZyYW1lKHlfYmMsIFgpKQogIHByZWRfYmMgPC0gcHJlZGljdChsbV9iYywgbmV3ZGF0YSA9IGRhdGEuZnJhbWUoWCkpCiAgcm1zZV9iYyA8LSBzcXJ0KG1lYW4oKHlfYmMgLSBwcmVkX2JjKV4yKSkKICAKICAjIFllby1Kb2huc29uIFRyYW5zZm9ybWF0aW9uCiAgeWpfb2JqIDwtIHllb2pvaG5zb24oeV9yYXcpCiAgeV95aiA8LSBwcmVkaWN0KHlqX29iaikKICAKICBsbV95aiA8LSBsbSh5X3lqIH4gLiwgZGF0YSA9IGRhdGEuZnJhbWUoeV95aiwgWCkpCiAgcHJlZF95aiA8LSBwcmVkaWN0KGxtX3lqLCBuZXdkYXRhID0gZGF0YS5mcmFtZShYKSkKICBybXNlX3lqIDwtIHNxcnQobWVhbigoeV95aiAtIHByZWRfeWopXjIpKQogIAogIHJlc3VsdHMgPC0gZGF0YS5mcmFtZSgKICAgIFRyYW5zZm9ybWF0aW9uID0gYygiTG9nIChsb2cxcCkiLCAiQm94LUNveCIsICJZZW8tSm9obnNvbiIpLAogICAgUk1TRSA9IGMocm1zZV9sb2csIHJtc2VfYmMsIHJtc2VfeWopLAogICAgQmVzdF9MYW1iZGFfQm94Q294ID0gYyhOQSwgYmVzdF9sYW1iZGEsIE5BKQogICkKICAKICBwcmludChyZXN1bHRzKQogIAogIHJldHVybihsaXN0KAogICAgbG9nX21vZGVsID0gbG1fbG9nLAogICAgYm94Y294X21vZGVsID0gbG1fYmMsCiAgICB5al9tb2RlbCA9IGxtX3lqLAogICAgYm94Y294X2xhbWJkYSA9IGJlc3RfbGFtYmRhLAogICAgdHJhbnNmb3JtYXRpb25fcmVzdWx0cyA9IHJlc3VsdHMKICApKQp9CmBgYAoKYGBge3J9CiMgU2VsZWN0IHRhcmdldCB2YXJpYWJsZQp0YXJnZXRfdmFyIDwtICJjcmFzaF9yYXRlX3Blcl8xMDAwIgoKZGYgPC0gZGYgJT4lIGRwbHlyOjpzZWxlY3QoLWluanVyeV9yYXRlX3Blcl8xMDAwLCAtZmF0YWxpdHlfcmF0ZV9wZXJfMTAwMCwgIC10b3RhbF9wb3B1bGF0aW9uLCAtZ2VvaWQpCgpmZWF0dXJlcyA8LSBzZXRkaWZmKG5hbWVzKGRmKSwgYygiY3Jhc2hfcmF0ZV9wZXJfMTAwMCIsICJib3JvdWdoIikpCmNvbXBhcmlzb24gPC0gY29tcGFyZV90cmFuc2Zvcm1hdGlvbnMoZGYsICJjcmFzaF9yYXRlX3Blcl8xMDAwIiwgZmVhdHVyZXMpCgpgYGAKCiMjIExPRyBUUkFOU0ZPUk0gVEFSR0VUCgpgYGB7ciBzaW5nbGUtdGFyZ2V0LWRmfQpYIDwtIGRmICU+JSBkcGx5cjo6c2VsZWN0KC10YXJnZXRfdmFyLCAtYm9yb3VnaCkKCiMgTG9nLXRyYW5zZm9ybSB0aGUgdGFyZ2V0IHZhcmlhYmxlCnkgPC0gbG9nMXAoZGZbW3RhcmdldF92YXJdXSkgIAoKaGlzdCh5LCBicmVha3MgPSA1MCwgbWFpbiA9ICJMb2ctVHJhbnNmb3JtZWQgQ3Jhc2ggUmF0ZSIsIHhsYWIgPSAibG9nMXAoY3Jhc2hfcmF0ZV9wZXJfMTAwMCkiKQpgYGAKCgojIyBIWVBFUlBBUkFNRVRFUiBUVU5JTkcgd2l0aCBPcHR1bmEKCmBgYHtyIHR1bmUtcGFyYW1zLCBldmFsPUZBTFNFfQojIyBDT05WRVJUIFRPIERNQVRSSVgKZHRyYWluX2FsbCA8LSB4Z2IuRE1hdHJpeChkYXRhID0gYXMubWF0cml4KFgpLCBsYWJlbCA9IHkpCgojIERlZmluZSBTcGF0aWFsIEZvbGRzIApib3JvdWdocyA8LSB1bmlxdWUoZGYkYm9yb3VnaCkKZm9sZHMgPC0gbGFwcGx5KGJvcm91Z2hzLCBmdW5jdGlvbihiKSB3aGljaChkZiRib3JvdWdoICE9IGIpKQpuYW1lcyhmb2xkcykgPC0gYm9yb3VnaHMKCiMjIFN0YXJ0IFB5dGhvbiB2ZW52CnJldGljdWxhdGU6OnVzZV92aXJ0dWFsZW52KCJyLXJldGljdWxhdGUiLCByZXF1aXJlZCA9IFRSVUUpCgojIyBPUFRVTkEtQkFTRUQgU1BBVElBTCBDVgpvcHR1bmEgPC0gaW1wb3J0KCJvcHR1bmEiKQoKIyBPYmplY3RpdmUgRnVuY3Rpb24gCm9iamVjdGl2ZSA8LSBmdW5jdGlvbih0cmlhbCkgewogIHBhcmFtcyA8LSBsaXN0KAogICAgYm9vc3RlciA9ICJnYnRyZWUiLAogICAgb2JqZWN0aXZlID0gInJlZzpzcXVhcmVkZXJyb3IiLAogICAgZXZhbF9tZXRyaWMgPSAicm1zZSIsCiAgICBldGEgPSB0cmlhbCRzdWdnZXN0X2Zsb2F0KCJldGEiLCAwLjAxLCAwLjMsIGxvZyA9IFRSVUUpLAogICAgbWF4X2RlcHRoID0gdHJpYWwkc3VnZ2VzdF9pbnQoIm1heF9kZXB0aCIsIDMsIDEyKSwKICAgIG1pbl9jaGlsZF93ZWlnaHQgPSB0cmlhbCRzdWdnZXN0X2ludCgibWluX2NoaWxkX3dlaWdodCIsIDEsIDEwKSwKICAgIHN1YnNhbXBsZSA9IHRyaWFsJHN1Z2dlc3RfZmxvYXQoInN1YnNhbXBsZSIsIDAuNSwgMS4wKSwKICAgIGNvbHNhbXBsZV9ieXRyZWUgPSB0cmlhbCRzdWdnZXN0X2Zsb2F0KCJjb2xzYW1wbGVfYnl0cmVlIiwgMC41LCAxLjApLAogICAgZ2FtbWEgPSB0cmlhbCRzdWdnZXN0X2Zsb2F0KCJnYW1tYSIsIDAsIDEwKSwKICAgIGxhbWJkYSA9IHRyaWFsJHN1Z2dlc3RfZmxvYXQoImxhbWJkYSIsIDAsIDEwKSwKICAgIGFscGhhID0gdHJpYWwkc3VnZ2VzdF9mbG9hdCgiYWxwaGEiLCAwLCAxMCkKICApCiAgCiAgY3ZfcmVzdWx0cyA8LSB4Z2IuY3YoCiAgICBwYXJhbXMgPSBwYXJhbXMsCiAgICBkYXRhID0gZHRyYWluX2FsbCwKICAgIG5yb3VuZHMgPSAxNTAsCiAgICBmb2xkcyA9IGZvbGRzLAogICAgZWFybHlfc3RvcHBpbmdfcm91bmRzID0gMjAsCiAgICB2ZXJib3NlID0gMCwKICAgIG50aHJlYWQgPSBkZXRlY3RDb3JlcygpIC0gMQogICkKICAKICByZXR1cm4oY3ZfcmVzdWx0cyRldmFsdWF0aW9uX2xvZyR0ZXN0X3Jtc2VfbWVhbltjdl9yZXN1bHRzJGJlc3RfaXRlcmF0aW9uXSkKfQoKCnNldC5zZWVkKDIwMjUpCnN0dWR5IDwtIG9wdHVuYSRjcmVhdGVfc3R1ZHkoZGlyZWN0aW9uID0gIm1pbmltaXplIikKc3R1ZHkkb3B0aW1pemUob2JqZWN0aXZlLCBuX3RyaWFscyA9IDUwKQoKYmVzdF9wYXJhbXMgPC0gc3R1ZHkkYmVzdF9wYXJhbXMKYmVzdF9wYXJhbXMkb2JqZWN0aXZlIDwtICJyZWc6c3F1YXJlZGVycm9yIgpwcmludChiZXN0X3BhcmFtcykKYGBgCiRldGEKWzFdIDAuMjU4NjI5NQoKJG1heF9kZXB0aApbMV0gMTIKCiRtaW5fY2hpbGRfd2VpZ2h0ClsxXSA0Cgokc3Vic2FtcGxlClsxXSAwLjk5MjE0MzUKCiRjb2xzYW1wbGVfYnl0cmVlClsxXSAwLjU4MjM4NTgKCiRnYW1tYQpbMV0gMC4wODMxNjkwNQoKJGxhbWJkYQpbMV0gNC4yNjA5MDcKCiRhbHBoYQpbMV0gMC4yMDQ5MTEyCgoKCiMjIFRSQUlOL1RFU1QgU1BMSVQKCmBgYHtyIHRyYWluLXRlc3Qtc3BsaXR9CiMgU2V0IHNlZWQKc2V0LnNlZWQoMjAyNSkKCiMgU3BsaXQgYnkgaW5kZXgKdHJhaW5faW5kZXggPC0gY3JlYXRlRGF0YVBhcnRpdGlvbih5LCBwID0gMC44LCBsaXN0ID0gRkFMU0UpClhfdHJhaW4gPC0gWFt0cmFpbl9pbmRleCwgXQp5X3RyYWluIDwtIHlbdHJhaW5faW5kZXhdClhfdGVzdCA8LSBYWy10cmFpbl9pbmRleCwgXQp5X3Rlc3QgPC0geVstdHJhaW5faW5kZXhdCgojIENvbnZlcnQgdG8geGdiLkRNYXRyaXgKZHRyYWluIDwtIHhnYi5ETWF0cml4KGRhdGEgPSBhcy5tYXRyaXgoWF90cmFpbiksIGxhYmVsID0geV90cmFpbikKZHRlc3QgPC0geGdiLkRNYXRyaXgoZGF0YSA9IGFzLm1hdHJpeChYX3Rlc3QpLCBsYWJlbCA9IHlfdGVzdCkKYGBgCgojIyBUUkFJTiBNT0RFTAoKYGBge3IgdHJhaW4tbW9kZWwsIGV2YWw9RkFMU0V9CiMgU2V0IHNlZWQKc2V0LnNlZWQoMjAyNSkKCiMgVHJhaW5pbmcgd2l0aCBwYXJhbGxlbCBwcm9jZXNzaW5nCmZpbmFsX21vZGVsIDwtIHhnYi50cmFpbigKICBwYXJhbXMgPSBsaXN0KAogICAgZXRhID0gYmVzdF9wYXJhbXMkZXRhLAogICAgbWF4X2RlcHRoID0gYmVzdF9wYXJhbXMkbWF4X2RlcHRoLAogICAgZ2FtbWEgPSBiZXN0X3BhcmFtcyRnYW1tYSwKICAgIGNvbHNhbXBsZV9ieXRyZWUgPSBiZXN0X3BhcmFtcyRjb2xzYW1wbGVfYnl0cmVlLAogICAgbWluX2NoaWxkX3dlaWdodCA9IGJlc3RfcGFyYW1zJG1pbl9jaGlsZF93ZWlnaHQsCiAgICBzdWJzYW1wbGUgPSBiZXN0X3BhcmFtcyRzdWJzYW1wbGUsCiAgICBvYmplY3RpdmUgPSAicmVnOnNxdWFyZWRlcnJvciIsCiAgICBldmFsX21ldHJpYyA9ICJybXNlIgogICksCiAgZGF0YSA9IGR0cmFpbiwKICBucm91bmRzID0gMTAwMCwKICB3YXRjaGxpc3QgPSBsaXN0KHRyYWluID0gZHRyYWluLCB0ZXN0ID0gZHRlc3QpLAogIGVhcmx5X3N0b3BwaW5nX3JvdW5kcyA9IDIwLAogIHZlcmJvc2UgPSAxLAogIG50aHJlYWQgPSBkZXRlY3RDb3JlcygpIC0gMQopCmBgYAoKIyMgU0FWRSBNT0RFTAoKYGBge3Igc2F2ZS1tb2RlbCwgZXZhbD1GQUxTRX0KIyBDcmVhdGUgZGlyZWN0b3J5IGlmIGl0IGRvZXNuJ3QgZXhpc3QKaWYgKCFkaXIuZXhpc3RzKCIuLi9kYXRhL21vZGVscyIpKSB7CiAgZGlyLmNyZWF0ZSgiLi4vZGF0YS9tb2RlbHMiLCByZWN1cnNpdmUgPSBUUlVFKQp9CgojIFNhdmUgdGhlIGZpbmFsIFhHQm9vc3QgbW9kZWwKc2F2ZVJEUyhmaW5hbF9tb2RlbCwgZmlsZSA9ICIuLi9kYXRhL21vZGVscy94Z2JfbW9kZWwucmRzIikKCiMgU2F2ZSB0aGUgYmVzdCBwYXJhbWV0ZXJzCnNhdmVSRFMoYmVzdF9wYXJhbXMsIGZpbGUgPSAiLi4vZGF0YS9tb2RlbHMveGdiX2Jlc3RfcGFyYW1zLnJkcyIpCgpjYXQoIk1vZGVsIGFuZCBwYXJhbWV0ZXJzIHNhdmVkIHRvIC4uL2RhdGEvbW9kZWxzLyIpCmBgYAoKYGBge3IgbG9hZC1tb2RlbH0KIyBMb2FkIHRoZSBmaW5hbCBYR0Jvb3N0IG1vZGVsCmZpbmFsX21vZGVsIDwtIHJlYWRSRFMoIi4uL2RhdGEvbW9kZWxzL3hnYl9tb2RlbC5yZHMiKQoKIyBMb2FkIHRoZSBiZXN0IHBhcmFtZXRlcnMKYmVzdF9wYXJhbXMgPC0gcmVhZFJEUygiLi4vZGF0YS9tb2RlbHMveGdiX2Jlc3RfcGFyYW1zLnJkcyIpCgpgYGAKCiMjIE1PREVMIEVWQUxVQVRJT04KCmBgYHtyIGV2YWwtbW9kZWwsIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0V9CnNldC5zZWVkKDIwMjUpCgojIFByZWRpY3Qgb24gdGVzdCBzZXQKcHJlZHMgPC0gcHJlZGljdChmaW5hbF9tb2RlbCwgYXMubWF0cml4KFhfdGVzdCkpCgojIE1ldHJpY3MgCnJtc2UgPC0gc3FydChtZWFuKCh5X3Rlc3QgLSBwcmVkcyleMikpCm1hZSA8LSBtZWFuKGFicyh5X3Rlc3QgLSBwcmVkcykpCnIyIDwtIDEgLSAoc3VtKCh5X3Rlc3QgLSBwcmVkcyleMikgLyBzdW0oKHlfdGVzdCAtIG1lYW4oeV90ZXN0KSleMikpCgpjYXQoIk1vZGVsIEV2YWx1YXRpb24gTWV0cmljczpcbiIpCmNhdCgiICBSTVNFOiIsIHJtc2UsICJcbiIpCmNhdCgiICBNQUUgOiIsIG1hZSwgIlxuIikKY2F0KCIgIFLCsiAgOiIsIHIyLCAiXG5cbiIpCgojIFJlc2lkdWFscyAKcmVzaWR1YWxzIDwtIHlfdGVzdCAtIHByZWRzCnJlc2lkdWFsX2RmIDwtIGRhdGEuZnJhbWUoCiAgYWN0dWFsID0geV90ZXN0LAogIHByZWRpY3RlZCA9IHByZWRzLAogIHJlc2lkdWFscyA9IHJlc2lkdWFscwopCgojIFByZWRpY3RlZCB2cyBBY3R1YWwKcDEgPC0gcmVzaWR1YWxfZGYgJT4lCiAgZ2dwbG90KGFlcyh4ID0gYWN0dWFsLCB5ID0gcHJlZGljdGVkKSkgKwogIGdlb21fcG9pbnQoYWxwaGEgPSAwLjUpICsKICBnZW9tX2FibGluZShzbG9wZSA9IDEsIGludGVyY2VwdCA9IDAsIGNvbG9yID0gInJlZCIpICsKICB0aGVtZV9taW5pbWFsKCkgKwogIGxhYnModGl0bGUgPSAiUHJlZGljdGVkIHZzIEFjdHVhbCBDcmFzaCBSYXRlcyIsCiAgICAgICB4ID0gIkFjdHVhbCIsCiAgICAgICB5ID0gIlByZWRpY3RlZCIpCgojIFJlc2lkdWFscyB2cyBQcmVkaWN0ZWQgCnAyIDwtIHJlc2lkdWFsX2RmICU+JQogIGdncGxvdChhZXMoeCA9IHByZWRpY3RlZCwgeSA9IHJlc2lkdWFscykpICsKICBnZW9tX3BvaW50KGFscGhhID0gMC41LCBjb2xvciA9ICJibHVlIikgKwogIGdlb21faGxpbmUoeWludGVyY2VwdCA9IDAsIGxpbmV0eXBlID0gImRhc2hlZCIsIGNvbG9yID0gInJlZCIpICsKICB0aGVtZV9taW5pbWFsKCkgKwogIGxhYnModGl0bGUgPSAiUmVzaWR1YWxzIHZzIFByZWRpY3RlZCIsCiAgICAgICB4ID0gIlByZWRpY3RlZCIsCiAgICAgICB5ID0gIlJlc2lkdWFscyIpCgojIFJlc2lkdWFsIERlbnNpdHkgCnAzIDwtIGdncGxvdChyZXNpZHVhbF9kZiwgYWVzKHggPSByZXNpZHVhbHMpKSArCiAgICBnZW9tX2hpc3RvZ3JhbSgKICAgICAgICBhZXMoeSA9IC4uZGVuc2l0eS4uKSwgIAogICAgICAgIGJpbndpZHRoID0gMC4wNSwgCiAgICAgICAgZmlsbCA9ICJzdGVlbGJsdWUiLCAKICAgICAgICBjb2xvciA9ICJ3aGl0ZSIsCiAgICAgICAgYWxwaGEgPSAwLjcKICAgICkgKwogICAgZ2VvbV9kZW5zaXR5KAogICAgICAgIGNvbG9yID0gInJlZCIsIAogICAgICAgIGxpbmV3aWR0aCA9IDEsIAogICAgICAgIGFkanVzdCA9IDEuNSAgCiAgICApICsKICAgIHRoZW1lX21pbmltYWwoKSArCiAgICBsYWJzKAogICAgICAgIHRpdGxlID0gIlJlc2lkdWFsIERpc3RyaWJ1dGlvbiIsIAogICAgICAgIHggPSAiUmVzaWR1YWxzIiwgCiAgICAgICAgeSA9ICJEZW5zaXR5IiAgIAogICAgKQoKCmdnc2F2ZSgiLi4vcmVwb3J0L3Bsb3RzL3ByZWRpY3RlZF92c19hY3R1YWxfdmFsdWVzX3Bsb3QucG5nIiwgcDEsIHdpZHRoID0gMTAsIGhlaWdodCA9IDYsIGRwaSA9IDMwMCkKZ2dzYXZlKCIuLi9yZXBvcnQvcGxvdHMvcmVzaXN1YWxzX3ZzX3ByZWRpY3RlZF92YWx1ZXNfcGxvdC5wbmciLCBwMiwgd2lkdGggPSAxMCwgaGVpZ2h0ID0gNiwgZHBpID0gMzAwKQpnZ3NhdmUoIi4uL3JlcG9ydC9wbG90cy9yZXNpZHVhbF9kZW5zaXR5X3Bsb3QucG5nIiwgcDMsIHdpZHRoID0gMTAsIGhlaWdodCA9IDYsIGRwaSA9IDMwMCkKCiMgR3JpZCBQbG90cwpncmlkX3Bsb3QgPC0gcDEgKyBwMiArIHAzICsgcGxvdF9sYXlvdXQobmNvbCA9IDMpCgpnZ3NhdmUoCiAgIi4uL3JlcG9ydC9wbG90cy9yZXNpZHVhbHNfZGlhZ25vc3RpY19ncmlkLnBuZyIsCiAgZ3JpZF9wbG90LAogIHdpZHRoID0gMTgsIAogIGhlaWdodCA9IDYsCiAgZHBpID0gMzAwCikKYGBgCk1vZGVsIEV2YWx1YXRpb24gTWV0cmljczoKICBSTVNFOiAwLjMzOTYxNDYgCiAgTUFFIDogMC4yNDIyNTM4IAogIFLCsiAgOiAwLjQyNTI4NSAKIyMgU0hBUCBFWFBMQUlOQUJJTElUWQoKYGBge3Igc2hhcH0KIyBDb21wdXRlIFNIQVAgdmFsdWVzCnNoYXBfdmFsdWVzIDwtIHNoYXAudmFsdWVzKHhnYl9tb2RlbCA9IGZpbmFsX21vZGVsLCBYX3RyYWluID0gYXMubWF0cml4KFhfdHJhaW4pKQpzaGFwX2xvbmcgPC0gc2hhcC5wcmVwKHNoYXBfY29udHJpYiA9IHNoYXBfdmFsdWVzJHNoYXBfc2NvcmUsIFhfdHJhaW4gPSBhcy5tYXRyaXgoWF90cmFpbikpCgojIFNIQVAgc3VtbWFyeSBwbG90CnByaW50KHNoYXAucGxvdC5zdW1tYXJ5KHNoYXBfbG9uZykpCgpzaGFwIDwtIHNoYXAucGxvdC5zdW1tYXJ5KHNoYXBfbG9uZykKZ2dzYXZlKCIuLi9yZXBvcnQvcGxvdHMvc2hhcF9zdW1tYXJ5X3Bsb3QucG5nIiwgc2hhcCwgd2lkdGggPSAxMCwgaGVpZ2h0ID0gNiwgZHBpID0gMzAwKQpgYGAKCiMjIElOU1BFQ1QgVFJFRSBFTlNFTUJMRQoKYGBge3IgaW5zcGVjdC10cmVlfQp4Z2IucGxvdC5tdWx0aS50cmVlcyhtb2RlbCA9IGZpbmFsX21vZGVsKQpgYGAKCmBgYHtyLCBldmFsPUZBTFNFfQojIFNIQVAgRGVwZW5kZW5jZSBhbmQgSW50ZXJhY3Rpb24gUGxvdHMKCiMgSW50ZXJhY3Rpb24gdmFsdWVzCnNoYXBfaW50ZXJhY3Rpb25fdmFsdWVzIDwtIHByZWRpY3QoCiAgZmluYWxfbW9kZWwsCiAgYXMubWF0cml4KFhfdHJhaW4pLAogIHByZWRpbnRlcmFjdGlvbiA9IFRSVUUKKQoKc2F2ZVJEUyhzaGFwX2ludGVyYWN0aW9uX3ZhbHVlcywgIi4uL2RhdGEvbW9kZWxzL3NoYXBfaW50ZXJhY3Rpb25zLnJkcyIpCmBgYApgYGB7cn0KIyBEZWZpbmUgdGhlIHNwZWNpZmljIHRvcCBmZWF0dXJlcyB3ZSB3YW50IHRvIGFuYWx5emUKc2VsZWN0ZWRfZmVhdHVyZXMgPC0gYygKICAicG9zdF9wYW5kZW1pYyIsCiAgInBjdF9hZ2VfNjVfcGx1cyIsCiAgInBjdF9ncmFkdWF0ZV9kZWdyZWUiLAogICJwY3RfaW5fbGFib3JfZm9yY2UiLAogICJtZWRpYW5faW5jb21lIiwKICAicGN0X2NvbW11dGVfbWVkaXVtIgopCgojIENyZWF0ZSBwcmV0dHkgbmFtZXMgZm9yIHRoZSBzZWxlY3RlZCBmZWF0dXJlcwpmZWF0dXJlX25hbWVzIDwtIGMoCiAgIlByZSB2cyBQb3N0IFBhbmRlbWljIiwKICAiJSBQb3B1bGF0aW9uIDY1KyBZZWFycyIsCiAgIiUgd2l0aCBHcmFkdWF0ZSBEZWdyZWUiLAogICIlIGluIExhYm9yIEZvcmNlIiwKICAiTWVkaWFuIEluY29tZSIsCiAgIiUgd2l0aCBNZWRpdW0gQ29tbXV0ZSIKKQoKIyBWZXJpZnkgYWxsIHNlbGVjdGVkIGZlYXR1cmVzIGV4aXN0IGluIFhfdHJhaW4KbWlzc2luZ19mZWF0dXJlcyA8LSBzZXRkaWZmKHNlbGVjdGVkX2ZlYXR1cmVzLCBjb2xuYW1lcyhYX3RyYWluKSkKaWYgKGxlbmd0aChtaXNzaW5nX2ZlYXR1cmVzKSA+IDApIHsKICBzdG9wKHBhc3RlKCJUaGUgZm9sbG93aW5nIGZlYXR1cmVzIGFyZSBtaXNzaW5nIGZyb20gWF90cmFpbjoiLCAKICAgICAgICAgICAgIHBhc3RlKG1pc3NpbmdfZmVhdHVyZXMsIGNvbGxhcHNlID0gIiwgIikpKQp9CgojIEdlbmVyYXRlIFBEUCBwbG90cyB3aXRoIGVycm9yIGhhbmRsaW5nCnBkcF9wbG90cyA8LSBtYXAyKHNlbGVjdGVkX2ZlYXR1cmVzLCBmZWF0dXJlX25hbWVzLCBmdW5jdGlvbihmLCBuYW1lKSB7CiAgdHJ5Q2F0Y2goewogICAgcGQgPC0gcGRwOjpwYXJ0aWFsKAogICAgICBmaW5hbF9tb2RlbCwgCiAgICAgIHByZWQudmFyID0gZiwgCiAgICAgIHRyYWluID0gYXMubWF0cml4KFhfdHJhaW4pLCAKICAgICAgZ3JpZC5yZXNvbHV0aW9uID0gMzAKICAgICkKICAgIAogICAgYXV0b3Bsb3QocGQpICsKICAgICAgZ2d0aXRsZShuYW1lKSArCiAgICAgIHRoZW1lX21pbmltYWwoKSArCiAgICAgIHRoZW1lKAogICAgICAgIGF4aXMudGl0bGUgPSBlbGVtZW50X2JsYW5rKCksCiAgICAgICAgYXhpcy50ZXh0ID0gZWxlbWVudF90ZXh0KHNpemUgPSA4KSwKICAgICAgICBwbG90LnRpdGxlID0gZWxlbWVudF90ZXh0KHNpemUgPSAxMCwgaGp1c3QgPSAwLjUpLAogICAgICAgIHBsb3QubWFyZ2luID0gdW5pdChjKDEsIDEsIDEsIDEpLCAibGluZXMiKQogICAgICApCiAgfSwgZXJyb3IgPSBmdW5jdGlvbihlKSB7CiAgICBtZXNzYWdlKHBhc3RlKCJFcnJvciBnZW5lcmF0aW5nIFBEUCBmb3IiLCBmLCAiOiIsIGUkbWVzc2FnZSkpCiAgICByZXR1cm4oTlVMTCkKICB9KQp9KQoKIyBSZW1vdmUgYW55IE5VTEwgZW50cmllcyBmcm9tIGZhaWxlZCBwbG90cwpwZHBfcGxvdHMgPC0gcHVycnI6OmNvbXBhY3QocGRwX3Bsb3RzKQoKIyBDcmVhdGUgZGlyZWN0b3J5IGlmIGl0IGRvZXNuJ3QgZXhpc3QKaWYgKCFkaXIuZXhpc3RzKCIuLi9yZXBvcnQvcGxvdHMvcGRwIikpIHsKICBkaXIuY3JlYXRlKCIuLi9yZXBvcnQvcGxvdHMvcGRwIiwgcmVjdXJzaXZlID0gVFJVRSkKfQoKIyBTYXZlIGluZGl2aWR1YWwgcGxvdHMKd2FsazIocGRwX3Bsb3RzLCBzZWxlY3RlZF9mZWF0dXJlcywgfnsKICBpZiAoIWlzLm51bGwoLngpKSB7CiAgICBnZ3NhdmUoCiAgICAgIGZpbGVuYW1lID0gcGFzdGUwKCIuLi9yZXBvcnQvcGxvdHMvcGRwL3BkcF8iLCAueSwgIi5wbmciKSwKICAgICAgcGxvdCA9IC54LAogICAgICB3aWR0aCA9IDUsCiAgICAgIGhlaWdodCA9IDQsCiAgICAgIGRwaSA9IDMwMAogICAgKQogIH0KfSkKCiMgQ3JlYXRlIGFuZCBzYXZlIGdyaWQgcGxvdAppZiAobGVuZ3RoKHBkcF9wbG90cykgPiAwKSB7CiAgZ3JpZF9wbG90IDwtIGdyaWRFeHRyYTo6Z3JpZC5hcnJhbmdlKAogICAgZ3JvYnMgPSBwZHBfcGxvdHMsIAogICAgbmNvbCA9IDIsCiAgICB0b3AgPSAiUGFydGlhbCBEZXBlbmRlbmNlIFBsb3RzIGZvciBLZXkgRmVhdHVyZXMiCiAgKQogIAogIGdnc2F2ZSgKICAgIGZpbGVuYW1lID0gIi4uL3JlcG9ydC9wbG90cy9wZHAvcGRwX2dyaWQucG5nIiwKICAgIHBsb3QgPSBncmlkX3Bsb3QsCiAgICB3aWR0aCA9IDEwLAogICAgaGVpZ2h0ID0gMTIsCiAgICBkcGkgPSAzMDAKICApCiAgCiAgIyBPcHRpb25hbGx5IGRpc3BsYXkgdGhlIGdyaWQgcGxvdAogIHByaW50KGdyaWRfcGxvdCkKfSBlbHNlIHsKICB3YXJuaW5nKCJObyBQRFAgcGxvdHMgd2VyZSBzdWNjZXNzZnVsbHkgZ2VuZXJhdGVkLiIpCn0KYGBgCgoK